% Julien Meaud
% revised by Thomas Bowling
%
% 11/13/2014
% compute the FTT of time domain simulations
% w should be the matrix of the BM velocity or displacement
%           1st index of w: cross-section
%           2nd index of w: time
% 
% 11/17/2017: Applies a flat top window
%-------------------------------------------------------------------------


function [vBMfreq,freq]=computeFFT2_FlatTopWindow(w,T,T0,Tf)
% vBMfreq(x,freq)

[N,N2]=size(w);
if(N > 1 && N2 > 1)
    
    deltaT=T(2)-T(1);
    
    % i0 = round(T0/deltaT) + 1;
    % i1 = round(Tf/deltaT) + 1;
    [i0,~,~] = find2(T, T0);
    [i1,~,~] = find2(T, Tf);
    
    
    
    if(T(1) ~= 0)
        i0 = 1;
        i1 = min(N2,length(T));
        
    end
    
    Ttotal=T(i1)-T(i0);  %total recorded time
    Nfft=i1-i0+1; %total number of samples
    fS=1/deltaT;  %sampling frequency
    deltaF=1/Ttotal;
    
    fN=fS/2;  %Nyquist frequency
    freq=0:deltaF:fN;

    
    
    FTwindow=flattopwin(Nfft,'periodic');
    FTwindow=FTwindow/(sum(FTwindow)/Nfft);
    
    vBMfreq = zeros(N,length(freq));
    
    %% Compute FFT for each position along BM
    for iX=1:N
        vBM=w(iX,i0:i1);
        if(~all(size(vBM)==size(FTwindow)))
            FTwindow=transpose(FTwindow);
        end
        
        V=2*fft(vBM.*FTwindow)/(double(Nfft));
        vBMfreq(iX,:)=V(1:length(freq));
    end
    
    %% Round freq to get correct frequency values
    dFFT = 1/(Tf-T0);
    freq = round(freq./dFFT) .* dFFT;
else
    [freq,vBMfreq] = computeFFT3_FlatTopWindow(w,T,T0,Tf);
end
end